/*
    Copyright 2010-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Implementations of camera projections.

// $Id$

#include "projection.h"
#include <CCfits/CCfits>
#include <exception>

using namespace std;
using CCfits::ExtHDU;

const string mcrx::rectilinear_projection::typestring_="rectilinear";
const string mcrx::aitoff_projection::typestring_="aitoff";

boost::shared_ptr<mcrx::projection> 
mcrx::projection::projection_factory(const string& type,
				     T_float fov, T_float aspect)
{
  if (type==rectilinear_projection::typestring_)
    return boost::shared_ptr<mcrx::projection>( new rectilinear_projection(fov, aspect));
  else if (type==aitoff_projection::typestring_)
    return boost::shared_ptr<mcrx::projection>( new aitoff_projection());

  cerr << "Unknown projection\n";
  throw runtime_error("unknown projection");
}

boost::shared_ptr<mcrx::projection> 
mcrx::projection::read_projection(ExtHDU& hdu)
{
  string camera_type; 
  hdu.readKey("camtype", camera_type);

  if (camera_type==rectilinear_projection::typestring_)
    return boost::shared_ptr<mcrx::projection>( new rectilinear_projection(hdu));
  else if (camera_type==aitoff_projection::typestring_)
    return boost::shared_ptr<mcrx::projection>( new aitoff_projection());

  cerr << "Unknown projection reading FITS file\n";
  throw runtime_error("unknown projection");
}


/** Projects a ray origin on a pixel in the image. Point is in a frame
    where the camera is pointing in the positive z-direction. */
blitz::TinyVector<mcrx::T_float,2>
mcrx::rectilinear_projection::project (vec3d point) const
  // point intentionally passed by value
{
  const bool abort_on_miss (false);
  const bool check_for_behind (true);

  // if point[2] < 0, we have to propagate the ray antiparallel to the
  // camera axis to get to the camera position, which means that the
  // ray position is behind the camera. Without this check, emission
  // behind the camera is also imaged, mirror-flipped, which is not
  // good. (In most "scientific" simulations, this will not matter,
  // but if we want to fly through the emitting volume it does.)
  if(check_for_behind && (point[2]<0)) {
    if (abort_on_miss)  
      assert (0); 
    return blitz::TinyVector<T_float,2> (blitz::quiet_NaN(T_float()),
					 blitz::quiet_NaN(T_float()));
  }
  
  /*
  // ignore emission coming from very close to the camera (only used
  // for movies)
  if(dot(point,point)<1) {
    return blitz::TinyVector<T_float,2> (blitz::quiet_NaN(T_float()),
					 blitz::quiet_NaN(T_float()));
  }
  */
  
  // Calculate "path length" (taking into account the non-normalized
  // dir) from ray position to the projection plane. (lpp-point).lpn is
  // the distance parallel to the camera axis, and lpn.dir is the
  // direction cosine along said axis.
  const T_float l=(point[2]-lpp_[2])/point[2];
  point *= 1-l;

  // (Note that l can be negative, that's normal because for very
  // small fields of view, the projection plane can actually be far in
  // front of the camera)
  
  // The ray is now in the projection plane. Calculate the components
  // of position in that plane.

  const T_float ix=point[0];
  const T_float iy=point[1];

  // Check if it lands within the image
  if( (abs(ix)> aspect_) || (abs(iy)> 1) ) {
    if (abort_on_miss)  
      assert (0); 
    return blitz::TinyVector<T_float, 2> (blitz::quiet_NaN(T_float()),
					 blitz::quiet_NaN(T_float()));
  }
  return blitz::TinyVector<T_float, 2> ((ix+aspect_)/(2*aspect_),(iy+1)*0.5);
}


/** Projects a ray origin on a pixel in the image. Point is in a frame
    where the camera is pointing in the positive z-direction located
    at the origin. Point is intentionally passed by value. */
blitz::TinyVector<mcrx::T_float,2>
mcrx::aitoff_projection::project (vec3d point) const
{
  const bool abort_on_miss (false);
  const bool check_for_behind (true);

  // Calculate the vector from ray position to camera (non-normalized
  // direction vector, but it doesn't matter), and rotate it into a
  // frame where the camera direction is the z axis and nx_ and ny_
  // are the x and y axes.

  // Calculate the "longitude" (angle in xz plane) and "latitude"
  // (angle in yz plane)
  const T_float lon = atan2(point[0],point[2]);
  const T_float lat = atan(point[1]/sqrt(point[0]*point[0]+point[2]*point[2]));

  // now map lat and long to (+-2^1.5,+-2^0.5) with to Hammer-Aitoff projection

  const T_float denom = sqrt(1+cos(lat)*cos(0.5*lon));

  const T_float ha_x = 2*sqrt(2.)*cos(lat)*sin(0.5*lon)/denom;
  const T_float ha_y = sqrt(2.)*sin(lat)/denom;

  // map to (0,1). (The fixed aspect will restore area-preserving property.)
  return blitz::TinyVector<T_float, 2> (ha_x/(4*sqrt(2.))+0.5,
					ha_y/(2*sqrt(2.))+0.5);
}


/** Creates a camera projection by reading parameters from HDU keywords. */
mcrx::rectilinear_projection::rectilinear_projection (CCfits::ExtHDU& hdu) :
  fov_(1), aspect_(1),
  lpn_ (vec3d(0,0,-1)),
  ny_(vec3d(0,1,0)),
  nx_ (cross (ny_, lpn_)) 
{
  // Because the quantities are const, we have to resort to this
  // ugly-ass thing which might not even be valid.
  hdu.readKey("fov", const_cast<T_float&>(fov_));
  hdu.readKey("aspect", const_cast<T_float&>(aspect_));
  const_cast<vec3d&>(lpp_) = vec3d(0,0,1./tan(0.5*fov_));
}

/** Writes camera projection parameters as keywords to a FITS HDU. */
void 
mcrx::rectilinear_projection::write_parameters (CCfits::ExtHDU& hdu,
						const T_unit_map& u) const
{
  // write all the quantities as keywords to the specified HDU
  hdu.addKey("camtype", type (), "Camera model used to project rays");
  hdu.addKey("fov", fov_,
	     "[rad] Field of view");
  hdu.addKey("aspect", aspect_, "[] Image aspect ratio");
  hdu.addKey("solid_angle", solid_angle(),
	     "[sr] Image area in steradians");
}



/** Writes camera projection parameters as keywords to a FITS HDU. */
void 
mcrx::aitoff_projection::write_parameters (CCfits::ExtHDU& hdu,
					   const T_unit_map& u) const
{
  // write all the quantities as keywords to the specified HDU
  hdu.addKey("camtype", type (), "Camera model used to project rays");

  hdu.addKey("solid_angle", solid_angle(),
	     "[sr] Image area in steradians");
}

